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This paper examines the distribution of the mantissas of floating point 
numbers and shoios hoiv the arithmetic operations of a computer transform 
various distributions toward the limiting distribution 

(where b is the base of the number system) . The paper also gives a number 
of applications to hardware, software, and general computing which show 
that this distribution is not merely an amusing curiosity. A brief exami- 
nation of the distribution of exponents is included. 

I. INTRODUCTION 

The main purpose of this paper is to examine, from the computing 
machine's point of view, the well-known (to comparatively few people) 
unequal distribution of the "mantissas" of "naturally occurring" sets 
of numbers. The observed probability density distributions are often 
close to the reciprocal density distribution 
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where b is the number base (usually 2, 8, 10, or 16). The corresponding 
cumulative probability distribution is 



R(t) = f r(x) civ = f -~- 
Ji/ b Ji/b x In 



where, of course, 



In t + In b 
In b 



12(1/6) = and 12(1) - 1. 



(2) 



From the cumulative distribution, it follows that the probability of 
observing the leading digit N of a number that is drawn at random 
from r(t) is 



R(N + 1) - R(N) = 



In (JV + 1) - In (N) 
In b 



(3) 



and this is usually what is measured in experiments. 

A typical experiment is that of tabulating the number of physical 
constants in a table having a given leading digit (see Table I and Ref. 1, 
p. 7). The result looks reasonable. Many other examples of observing 
the reciprocal distribution have been reported. For references see 
Refs. 2 and 3. 

The reciprocal distribution has been explained in many ways. One 
popular but not immediately obvious explanation for the distribution of 
physical constants is as follows. Consider the distribution of the leading 



Table I — The Distribution of the Leading 
Digits of 50 Physical Constants 





Number of cases 


Expected number 




Leading digit N 


observed 


eq. (3) 


Difference 
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16 


15 
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11 


5) 
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-4 
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5 


5 
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2 
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-1 
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3 


—2 
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3 


2 


1 




50 


50 
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digits of the set of all the physical constants that might occur. If the 
units of measurement were to be changed then the corresponding leading 
digit of any particular physical constant would probably change, but 
it is difficult to believe that the distribution itself would change sig- 
nificantly. To believe so seems to indicate a belief that either the present 
units of measurement or else the new set have some intimate connection 
with the real world. An alternative, and more elegant, explanation is 
given by Roger Pinkham in his classic paper (Ref. 2). The explanation 
given in the present paper is based on how the computer transforms 
distributions during arithmetic operations. In particular the paper shows 
how, from any reasonable distributions, repeated multiplications and/or 
divisions rapidly move the distributions toward the reciprocal dis- 
tribution. The effect for addition and subtraction is somewhat different. 
The paper also shows the persistance of the reciprocal distribution 
once it is attained. 

Since floating point numbers are the basis of most of numerical 
analysis one may well ask why this obvious and experimentally well- 
verified distribution is so often ignored. Is it because it appears to 
contradict the usually accepted model of the number system in which 
numbers correspond to points on a homogeneous straight line? Not only 
are the floating point numbers not uniformly spaced in a computer 
(the difference between the two largest possible numbers is very large, 
while the distance between the two smallest positive number is very 
small, and zero is relatively isolated), but the reciprocal distribution 
shows that even in intervals in which the numbers are equally spaced 
they are not equally likely to occur. 

Thus in analogy with non-Euclidean geometry this paper proposes 
an alternative to the conventional identification of numbers with points 
on a homogeneous straight line. Instead of adopting a measure for sets 
that is invariant under translation 

x' = x + k, 

we often prefer a measure that is invariant under scaling, namely 

x' = kx (k ^ 0). 

The reciprocal distribution is of practical as well as theoretical interest 
as we shall show in Section VII. In view of these examples, it is hoped 
that by adopting the machine's point of view with respect to how 
numbers are transformed by arithmetical operations, the computer 
scientists will become more aware of the importance of this distribution 
in manj' situations including numerical analysis. 
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II. THE MODEL 

The floating point numbers in a computing machine form a discrete, 
finite set. As is true in so many applications of mathematics to practical 
problems, we shall approximate a discrete distribution by a continuous 
one of sufficient smoothness. Anyone familiar with the upper and lower 
Riemann Integral sums can appreciate the degree of approximation 
being made (provided common sense is used in choosing the values of 
the curve between the given points). In the limit of the Riemann sum 
all the | Ax t | become less than any given e > 0; we of course need to 
stop at the granularity of the number system used, typically 10~ 8 or 
smaller. 

In principle, it is possible to carry this error estimate throughout all 
the subsequent steps of the mathematics to see how much the mathe- 
matics errs from reality; but it is customary to recognize that a little 
intuition will suffice to convince the user that the error will be much less 
than the accuracy of the experiments that the theory is designed to 
account for. Thus we have no need to get excited about such things as 
the Banach measure of a set (Ref . 4) ; we do not intend in this paper to 
let the mathematics obscure what is going on. The fact that computers 
are finite and operate at a finite speed for a finite length of time spares us 
from taking seriously all the confusions that can arise in mathematics 
when dealing with the infinite. 

III. THE BASIC FORMULAS 

In this section we derive the basic formulas which describe how dis- 
tributions are combined and transformed by the four arithmetic opera- 
tions of a computer. Let f(x) be the density distribution of the factor x, 
g{y) be the density distribution of the factor y, and h(z) be the density 
distribution of the result z of the arithmetic operation. Further, let 
F(x), G(y), and H(z) be the corresponding cumulative distributions. 

For both multiplication and division, the mantissas are directly 
combined and the exponents do not enter into the formation of the 
distribution of the result of the operation. Thus, it is sufficient in these 
cases to consider the distributions for (1/6 ^ x, y ^ 1). 

For multiplication, an examination of Fig. 1 shows that when the 
product falls in the shaded regions then the mantissa of the product is 
in the interval (1/fr, z). Thus the cumulative distribution H(z) is given by 

H(z) =['['' f(x)g(y) dy dx + f f f(x)g(y) dy dx 

•>l/b J\/b J\/b Jl/bx 

+ / / Kx)g(y) dy dx 

Jm Jl/bx 
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= f 1(x)[G(z/bx) - G(l/b) + G(l) - G(l/bx)] dx 

•'1/6 

+ £ Kx)[G(z/x) - G(l/bx)] dx. 

Differentiating with respect to z to get the density distribution we have 
h(z) = 1(z)[G(l/b) - (7(1/6) + (?(1) - (7(1/68) - (?(1) + G(l/bz)] 

+ /"' 1(x)g(z/bx){l/bx) dx + f Kz)g<j./z)Q./z) dx 

= \f M g(*/bx) dx + f ^ ^/x) <fe. (4) 

Similarly for division. The shaded region of Fig. 2 shows where the 
quotient x/y is less than z; thus the cumulative distribution for the 
quotient is 

H(z) = f f /(x)fl(y) dy da; + f f Kx)g(y) dy dx 

J\/b J\/b J\/b ->z/z 



+ / / f(x)g(y)dydx 

Jz Jx/bz 



= f /(z)[(7(z) - (7(1/6) + (7(1) - (?(*/«)] dx 

•>l/b 

+ J f(x)[G(x) - G(x/bz)] dx. 
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Fig. 1 — The cumulative probability distribution for the product z = xy. 
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y=v y = x 




Fig. 2 — -The cumulative probability distribution for the quotient z = x/y. 

Again differentiating with respect to z to get the density distribution 
we have 

h(z) = f(z)[G(z) - G(l/b) + (7(1) - 0(1) - G{z) + (7(1/6)] 

+ f /(*)[-0(*/«)(-*/^]<k+ f K*)l-g(x/l*){-xW\dx 

J l/b •> z 

= \ f xf(x)g(x/z) dx + jjp £ xf(x)g(x/bz) dx. (5) 

For both addition and subtraction the difference in the exponents of 
the two numbers x and ?/ is used to shift one mantissa with respect to 
the other before they are combined. For addition, we may suppose that 
one of the numbers, say x, lies in the range z/2 ^ x S z. The other 
term, y, therefore lies in the range z/2 ^ y ^ z ■&"*, where k is the 
number of digits in the mantissa and we set b~ k = e. Thus the density 
distribution of the sum is 



h(z) = 



j{x)g{z - x) dx. 



(0) 



For subtraction we suppose, without loss of generality, that x ^ y > 0, 
and 

z = x — y 
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with z ^ x ^ z/t. Then the density distribution is given by 

h(z) = f x)g(z + x) dx. (7) 

We have now derived the basic relations for the density distributions 
that arise from combining two numbers from arbitrary distributions 
according to the four arithmetic operations of a computer. 

IV. THE PERSISTENCE OF THE RECIPROCAL DISTRIBUTION" 

111 this section, we first show for both multiplication and division 
that if one of the factors x or y comes from the reciprocal distribution, 
and regardless of the distribution of the other factor, then h(z) is the 
reciprocal distribution. In particular, if a number is chosen from the 
reciprocal distribution, then its reciprocal is also from the reciprocal 
distribution. For addition and subtraction we show somewhat less. 

For the product set 

f® - ih »> 

in equation (4). We get for any distribution f(x) 

... If /(..•) bx . f x 1(x) x 

h(z) = t I ■ : — r dx + / : — j dx 

b J 1/b x z In b J : x z In b 

l 



z In b 



/(.»•) dx + /(.r) dx 



z In b 



(9) 



Obviously since z = xy, the same applies if we assume that f(x) is the 
reciprocal distribution. 

For the quotient, again assume equation (S) and put it in equation (5). 



h ® - ? L ■ r/(J ' ) ^rti dx + b I *& dh * 
- rti {.L Kx) " r + f. Kx) *} = ih- 



(10) 



In the special case of j(x) being the "spike distribution" with all of its 
probability at x = 1 we see that the reciprocal of a variable having the 
reciprocal distribution has the reciprocal distribution. The case of x 
having the reciprocal distribution and producing the reciprocal distri- 
bution, regardless of the distribution of the denominator, is covered by 
the product form, or can be worked out directly if desired. 

Thus, if in a long sequence of multiplications and divisions at least one 
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factor has the reciprocal distribution, then regardless of how the distri- 
butions of the other factors are chosen the result is still the reciprocal 
distribution; the reciprocal distribution persists under multiplication 
and division and cannot be broken by any choices for the other factors. 
For addition let x come from the reciprocal distribution for some range 
with normalization factor N t , and y also come from a reciprocal dis- 
tribution with its corresponding range and normalization factor iV 2 . 
Then writing e = b~ k 



h(z) / — • 

Jz/1 - 1 " 



dx 



z — X 
f .u-«j x r x 

= N * N * ; ; + - 

■.'2/2 3 1 a; 



z 



fc] 





1 




T 

Z 


- 


.r_ 


|«(l-€ 






l«/a 







(11) 



where N 3 is some constant. 

Similarly for subtraction (different N t ) 



/e N, N 2 



»»-r? 



dx 



. SiN , V" 1 ("I _ 1 " 

J, z \x z + .r. 

M I'" 



da: 



^^ln 



= ^. (12) 

z 

It should be noted, however, that in the last two cases the assumption 
of the reciprocal distribution for such great ranges is suspicious to say 
the least, since we know from experience that all exponents are not 
equally likely. That the reciprocal distribution over a large range implies 
the equally likely distribution of the relevant exponents can be seen by 
examining the base 16 number system in exponents, but where the 
mantissas are in binary. Thus the mantissas can have one of the forms: 

O.lzxz- • • 

O.Olzx- • • 

O.OOlx--- 

0.0001- •• . 
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If we assume 



p(x) = 



1 



(-h^x^ 1), 



x In 16 

what are the probabilities of each of the four forms? For the first one 
1 . 1 



/ 



dx = 



j x In 16 4 In 2 



[ml 



Similarly, each of the others is \. This result is quite different from that 
of the flat distribution (see Table II). 

V. THE APPROACH TO THE RECIPROCAL DISTRIBUTION 

Having shown that once it arises the reciprocal distribution persists 
for multiplication and division, we need to show how it can arise. For 
this we need a measure of how far a distribution h(z) is from the recip- 
rocal distribution r(z). It is obvious that 



f [h(z) - r(z)} dz = 



(13) 



for any h(z) and this does not provide a useful measure of distance. We 
shall define the distance of h(z) from the reciprocal distribution r(z) by 



max 



Hz) - r(z) 
r(z) 



m D{h(z)} = D[h), 



(14) 



which measures the maximum of the difference relative to the reciprocal 
distribution (it is natural to use the relative error when dealing with 
floating point numbers). 



Table II — Probability of Observing Mantissas With Leading 
Zeros in Base 16 Numbers When Written in Base 2 





Range 


Binary 
Exponent 


Probabilities 


Form 


Flat 


Reciprocal 


0.0001 

O.OOlx 

O.Olxx 

O.lxxx 


1/16 g x g 1/8 
1/8 ^ x ^ 1/4 
1/4 <, x g 1/2 
1/2 £ z £ 1 


-3 

_2 

-1 




1/15 

2/15 
4/15 
8/15 


1/4 
1/4 
1/4 
1/4 
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We showed in equation (9) that for a product, 

r(z) = \ f %&r(z/bx) dx + [ ^r(z/x) dx. 

J\/b X J z X 

Subtracting this from equation (4) and dividing by r(z) we have 

Hz) - r(z) = 1 C' m \ g(z/bx) - r(z/bx) l & 
r(z) b J l/b x L r(z) J 

+ rm [ >/*) -rwo ] dx . 

J, x L r(z) J 



But 



bx 
bxr(z) = — , -' , = r(z/bx) 
z In o 



»* ) = *~if"& = r(V:r) ' 



and we have 



Hz) 



A/» L r(2/6a:) J 



+ M lfe *te? a >- (15) 

Since /(a;) ^ for (1/6 git^l), 



W ~ r(z) 
r(z) 



^ f 1(x)D{g)dx+ f j(x)D{g}clx 



^ Dig} 
for all z. From this it follows that 

D{h] ^ D{g\ (16) 

regardless of the choice of j{x). 

We note that the equality would hold if f(x) were a single spike at 
x = 1, say, but that in view of equation (13), we generally expect a 
great deal of cancellation in the square brackets of equation (15) as it 
is integrated over the range. 

It is easy to examine the rapidity of the approach in the case of all the 
factors coming from the flat distribution 

r \ 1 b 

p(x) = 



1-1/6 6-1 
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Equation (14) gives for two factors 

*> = HrH)Xf+fe)7: 

In the base b = 10, this is 

h(z) = £?{lnl0 - 9 1n«}, 



X 



(17) 



which (for the proper range) is given by Ref. 5 (p. 37). The distance of 
the flat distribution is 



max 

i/iogi§i 



10 In 10 
z-j- -1 



10 In 10 
9 



- 1 = 1.558... 



while the distance of equation (17) is equal to 0.3454- • • . See Table III 
for further results. 

Similarly for division using equations (10) and (5), we have 

■ g(x/z) -r(x/z) l ^ 
r(z) J 



Hz) - r(z) _. 1 



r{z) 



Z J 1/b 



But 







X 


= r(x/z) 




bz* r -& 
x 


= r(x/bz), 


and we have 








Hz) - r(z) 
r(z) 


£D[g)\ 


f Kx)da 

J\/b 



Table III — The Distance of a Continued Product as a 
Function of the Number of Factors Selected 
from a Flat Distribution 



Number of Factors 


Distance 


1 
2 
3 
4 


1.558 
0.3454 
0.0980 
0.0289 
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or 



D{h} ^D{g}. 
In the case of flat distributions 



h(z) = 



hjt+a 



2(6 - 
which for the base 10 is (see Ref. 5, p. 37) 

" (2) - s [ 10 + ?] 

and has a distance of 0.4071 • • • . 

For addition we select g(y) as a reciprocal distribution (with suitable 
normalization factor N), subtract the corresponding equations and 
divide by r(z) to get 



N* 



r< = ll L KX) *r 



Hz) - r(z) 



Ni_N 
X z 



-J[ 



/2 



«-* 



J AT 2 r(s)1 
|_z — z r(z)J 



dx 
r(z) 



r(x) 
But by the mean value theorem for integrals 

_ j z/2 z — x r{z) 



dx. 



Hz) - r{z) 
r(z) 



1(e)- 



_ tie) 



where z/2 ^ 6 ^ z(l — e). The integral has been shown in equation (11) 
to be exactly 1. Hence 

D{h(z)} £Dlf(x)}. 

A similar derivation works for subtraction. 

In view of the dubious assumption of having the reciprocal distribu- 
tion over a very large range we need to examine more carefully the 
behavior of the mantissas of sums of numbers selected from some 
distribution. Let us imagine a Monte Carlo experiment. We select 
numbers from the range (0 < a ^ x ^ b) having the probability density 
distribution p(x) with mean n and variance a 2 . Divide the range into 
n equal intervals 

(a, a + h), (a + h, a + 2h), ■ ■ • , [a + (n — l)h, b], 
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where h = (b — a)/n. By counting how many numbers fall in each 
interval we get estimates of p(x). 

Let us add 2'" numbers of this set of numbers. The range for the sum is 

(2 m a, 2 m 6), 

the mean mi = 2"Vi and a\ = 2V. But the central limit theorem says 
that the distribution of the sum approaches a normal distribution about 
the mean with half width a Y . Suppose, for convenience, that m fell in 
the middle of an interval. Then as m increases and we count the number 
of cases of mantissas in each interval (note that the m in the term 2 m 
appears in the exponent only) we will find more and more of them will 
fall in the interval containing y. (which has the same mantissa as /ti); 
the distribution approaches a spike! This does not contradict the central 
limit theorem; it merely says that if y. ^ (ji = is the exceptional case), 
the distribution contracts as seen from the point of view of floating 
point numbers. In loose words, standing at the origin and viewing the 
rapidly receding mean p, , the width of the distribution o-j seems to get 
narrower as compared to the sum — the sum recedes as 2 m , the half 
width changes as 2 m/2 . 

VI. RANGE OF EXPONENTS 

It is now clear that in order to examine carefully the effect of addition 
(and subtraction) on the reciprocal distribution, it is necessary to know 
the distribution of the exponents of the numbers to be combined. 
Unfortunately at this time about the only model we have is as follows. 
Assume a distribution of exponents. Under multiplication and division 
the exponents are added and subtracted (with, due to carries an extra 1 
occasionally added, or subtracted) and by the central limit theorem 
we can expect: (i) that the distribution of the exponents will approach 
a normal distribution (assuming that overflow and underflow do not 
happen first) and (it) that this distribution will gradually spread out 
proportional to the square root of the number of operations. Thus, it 
appears that in practice the distribution of exponents is probably not 
stationary. Addition tends to eliminate the smaller exponents, while 
subtraction tends to increase them. 

Experience in numerical analysis shows that the range of the output 
numbers is usually much greater than the range of the input numbers, 
enough so to make one suspect that the variance increases as indicated 
in the above model. 

As one thinks carefully about the matter of addition and subtraction 
it seems reasonable to believe that they will not greatly perturb the 
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reciprocal distribution; and the experimental data from "naturally 
occurring numbers", which must have included some additions and 
subtractions, seem to bear out this belief. 

The feeling that under repeated additions and subtractions the central 
limit theorem applies to numbers (which is true), and therefore con- 
tradicts the reciprocal distribution of the mantissas, is typical of the 
"fixed point arithmetic" viewpoint of numbers — we are representing the 
sums and differences as floating point numbers, and it is the distribution 
of these mantissas and their possible approach to the reciprocal distri- 
bution that is of relevance here. 

VII. APPLICATIONS OF THE RECIPROCAL DISTRIBUTION 

Besides accounting for the experimentally found distributions, the 
reciprocal distribution is relevant to many optimization situations. 

As a first example, 6 consider the problem of placing the decimal 
(binary) point in the number representation system in order to minimize 
the number of normalization shifts after the computation of a product. 
(It was probably the minimization of normalizing shifts that caused 
IBM to adopt the base 16 in the system 360). If the point is placed 
before the first digit, then products of the form 

O.xxx... 
O.xxx... 

O.Oxz... 

will require a shift to normalize the result; while if it is placed after the 
first digit, then products like 

x.x... 
x.x... 

xx.x... 

will require a shift. Clearly these two cases have complementary prob- 
abilities. For the reciprocal distributions the probability p of 

xy ^ 1/6 



is 



r 1 r l/bx _i l 

Ji/b Ji/b xlnbyh 



dy dx 



m r i (m i/fa - to iM !lx _ f i (ju) dx 

J l/b In b \ x I J i/b In b \ x I 



In 2 6 



In 2 * 
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But for a flat distribution, 
b 



V 



if L C dy dx = (b=i) L (i - v dx 



b - 

th)i\* h -V-i 



b In 6 - (6 - 1) 
(6 - D 2 



For 6 = 2 this is 



p = 2 In 2 - 1 ^ 0.38. 



As a second application, consider the estimation of the effect of the 
representation error of numbers in base 2 and base 16. In Ref. 7 
McKeeman reports that the maximum relative representation error 
(MRRE) and the average relative representation error (ARRE) are 
as shown in Table IV, where the average is over the reciprocal distri- 
bution. 

A third example is the application to roundoff propagation. If x x has 
an error t x and x 2 has error e 2 , then in the product 

Xi + €i 

xz -V to 



X&% + X& + X 2 € 1 + 6l€ 2 

it is the leading digits that control the estimate of the propagated error. 
For the reciprocal distribution the mean is 

1-1/6 6-1 



J i/b x In 6 



dx = 



In 6 



6 In 6 



For base 2, this is 



x = 



2 In 2 



S 0.72134. 



Table IV — Maximum Relative Representation Error and 
Average Relative Representation Error 





MRRE 


ARRE 


binary 


1/2 X 2"" 


0.18 X 2-" 


octal 


•J-37 


0.21 X 10-" 


hexadecimal 


2~37 


0.17 X 2-" 
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The second moment about the mean is 



If 

In b J u 



which for 6 = 2 is 



6 2 In 6 l 2 In b 



M ' - ih d - ss) = °- 020674 - 

For the flat distribution, x = 0.75 and Af 2 = 0.020833. 

Thus we see that the effect of the reciprocal distribution on the average 
roundoff propagation is surprisingly small. 

Another example in which the reciprocal distribution must be con- 
sidered is that of producing "random" floating point mantissas. To 
generate these mantissas we use the earlier result that a long sequence 
of multiplications of numbers from a flat distribution will approximate 
a reciprocal distribution. Thus random mantissas can be generated by 

Y n = Y n . l -r n (shifted) 

where r n is from the usual (fiat) random number generator and "shifted" 
means after each product the leading zeros are shifted off. How well does 
this work? Experimental verification* is given by 8192 trials. Counting 
the number of mantissas falling in each of N catagories (see Table V). 

The last two columns of Table V give the sign changes observed in 
the difference between the observed and theoretical reciprocal distri- 
bution. The expected number of sign changes might be expected to be 
(N — l)/2, but since for N = 2 it is clear that one sign change will 
occur (because the mean of the residuals is zero) we have used N/2 as 
the expected number. The chi-square test shows that the two distribu- 
tions are close; the sign change test shows that the residuals are not 
systematically distributed. From these tests, we see that the generator 
"works." It is interesting to note that the period of this generator may 
well be much longer than that of the underlying flat random number 
generator. 

It is easy to see as a general rule that when we try to optimize a 
library routine for minimum mean running time (as against the 
Chebyshev minimax run time) we need to consider the distribution of 
the input data. Hence floating point numerical routines need to consider 
the reciprocal distribution; the square root, log, exponential, and sine 



Thanks to Brian Kernighau. 
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Table V — Distribution of 8192 Random Mantissas 





X 2 


Degrees of 
Freedom 


Residuals 


N 


Sign 
Changes 


Expected 


64 


61.392 


63 


30 


32 


32 


22.804 


31 


14 


16 


16 


11.150 


15 


8 


8 


8 


7.724 


7 


5 


4 


4 


3.261 


3 


2 


2 


2 


1.467 


1 


1 


1 



are all examples. In the case of the exponential and sine, some study of 
the exponents is also necessary. 
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